arXiv:1507.08405vl [stat.ME] 30Jul2015 


Statistical Science 

2015, Vol. 30, No. 2, 181-183 

DOI: 10.1214/15-STS519 

Main article DOI: 10.1214/14-STS487 

© Institute of Mathematical Statistics, 2015 

Rejoinder 

Marc G. Genton and William Kleiber 


We are grateful to the discussants for providing 
very valuable and insightful comments. Next, we 
present our views on some of the comments of the 
discussants and provide further discussion. 

We thank Bevilacqua, Hering and Porcu (here¬ 
after, BHP) for bringing attention to the fundamen¬ 
tal problem of comparing multivariate models. Until 
now, almost all comparisons between models have 
been relegated to empirical performance on specific 
datasets, whether it be performance on cokriging 
or particular scoring rules. BHP introduce two the¬ 
oretical approaches to comparing the flexibility of 
multivariate frameworks: (A) assessing the size of 
allowable co-located cross-correlation between pro¬ 
cesses, and (B) a measure of difference in allowed 
spatial (cross)-correlation at differing distances. 

Regarding (A), BHP claim the bivariate Matern is 
less flexible than the LMC in that there are nontriv¬ 
ial restrictions on the cross-correlation coefficient for 
the bivariate Matern that are not present for the 
LMC. We emphasize, however, that the bivariate 
Matern restrictions are a characterizing feature of 
the covariance class—no LMC construction can al¬ 
low for marginal and cross exact Matern behavior 
while allowing for unrestricted choice of co-located 
cross-correlation. Rather, it is a physical restriction 
on the covariance class, not a flexibility restriction. 

Most spatial modelers include a nugget effect in 
the statistical model, K)(s) = Z*( s) + £j(s), where 
Zi( s) is endowed with a multivariate model, and 
£j(s) is a white noise process that is uncorrelated 
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with Zi( s). If £j(s) is nontrivial with variance r 2 , 
then the restrictions on the cross-correlation coef¬ 
ficient can be relaxed, the amount depending on 
the magnitude of the nugget effect and sample size. 
To see this, let p = 2 and write the covariance ma¬ 
trix for two unit variance processes at n locations 
{Zi(si),...,Zi(s n ),Z 2 (si),...,Z 2 (s n )} T as B © £, 
where £ = {Cij(s k ,s e )}^ =1 . k/=1 and B = (B i:j ) 2 l j=x 
consists of four n x n block matrices. For simplic¬ 
ity, assume t\ =t 2 , so that Ri 2 = B 2 1 are matrices 
populated by a constant po and Bn = R 22 are ma¬ 
trices of ones with diagonal 1 + r 2 . Note that the 
case po = 1 results in B0S having the specified 
multivariate dependence; if po > 1, then the two pro¬ 
cesses can have larger cross-correlation than allowed 
by the specified model. The cases where po > 1 are 
valid when B i 2 = B 1X KB 22 , where K is a contrac¬ 
tion matrix (i.e., a matrix whose singular values are 
bounded by unity); this follows from Proposition 1 
of Kleiber and Genton (2013). This is one feasible 
way to relax the restrictions that are suggested by 
BHP’s (A) criterion. We view BHP’s (B) as an al¬ 
ternative interesting route to comparing models, al¬ 
though it is still unclear what improvements a mod¬ 
eler would expect to gain for various magnitudes of 
the (B) criterion. 

Cressie et al. focus on three main aspects: the im¬ 
portance of modeling the nugget effect (which yields 
additional potential difficulties in the multivariate 
context), the pseudo cross-variogram and alterna¬ 
tive approaches to building multivariate structures. 

We focused our efforts on reviewing multivari¬ 
ate covariance functions, not multivariate model¬ 
ing, a byproduct of which is that we left little dis¬ 
cussion to the issue of modeling the nugget effect. 
For instance, the underlying latent smooth process 
W of Cressie et al. [(2015), equation (4)] still re¬ 
quires specification of the multivariate structure, re¬ 
gardless of whether a nugget effect will or will not 
ultimately be included. Nonetheless, these authors 
bring up an important point in that, especially for 
multivariate processes, some variables may be mea¬ 
sured by the same instrument, in which case it may 
be expected that measurement errors are correlated 
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across variables at individual locations. Disentan¬ 
gling microscale variability of the process from mea¬ 
surement error is indeed a difficult prospect; Sang, 
Jun and Huang (2011) used a full-scale approxima¬ 
tion for multivariate processes that explicitly breaks 
up large scale, small scale and measurement error 
variability. 

Cressie et al. champion a traditional geosta- 
tistical approach to estimation, using a weighted 
least squares distance from empirical pseudo cross- 
variograms to estimate parameters of a parametric 
class of cross-covariance functions. Although this 
is a feasible route to estimation, interpretation of 
the pseudo cross-variogram remains unclear, unless 
the process is stationary and the variables are stan¬ 
dardized, so that the pseudo cross-variogram can be 
rewritten Var{Zi(s + h)/ci — Z 2 (s)/cj 2 }, where cq 
is the standard deviation of the 1th process. This 
seems a particularly important consideration if the 
scales of the two processes differ by orders of magni¬ 
tude. A more modern approach to estimation (albeit 
one that requires more modeling assumptions) is 
to adopt a maximum likelihood or Bayesian frame¬ 
work. Even if the processes are not strictly Gaussian, 
say, these other approaches can still be used for esti¬ 
mation if the deviation from Gaussianity is not too 
great. 

We agree with Cressie et al. that the conditional 
approach to estimation has utility in certain situ¬ 
ations when there is clear directional dependence 
between variables. However, if many (say, greater 
than four) processes are considered simultaneously, 
we echo Cressie et al.’s caution: “when more vari¬ 
ables are involved, the order may not always be ob¬ 
vious but, if the goal is to construct valid covariance 
and cross-covariance functions, the different order¬ 
ings can be viewed as enlarging the space of valid 
models.” On the other hand, the factor process ap¬ 
proach seems promising. Given the difficulties basis 
decomposition approaches experience in the univari¬ 
ate setting (Finley et al. (2009), Stein (2014)), we 
expect similar issues to arise in the multivariate set¬ 
ting. Thus, we suggest a multiresolution approach, 
decomposing differing scales of support into various 
resolutions of basis functions (Nychka et ah, 2002, 
2015). 

Finally, Cressie et ah bring up some important is¬ 
sues in validating and comparing statistical models 
that we view as sensible guidelines for future authors 
working in this field. Indeed, direct likelihood com¬ 
parisons may be muddled by the varying numbers 


of parameters between models, and BHP have intro¬ 
duced some tools apart from the usual information 
criteria to compare multivariate models. 

We are pleased that Simpson, Lindgren and Rue 
(hereafter, SLR) decided to expand on our very brief 
mention of the spectral representation of the cross¬ 
covariance matrix function, which for simplicity was 
restricted to the symmetric case in our review. Al¬ 
though we did address the topic of asymmetric cross¬ 
covariance functions in Section 5.1, SLR now pro¬ 
vided information about the spectral representa¬ 
tion in the asymmetric setting and further discussed 
the spectral representation of multivariate Gaus¬ 
sian random fields themselves. SLR also argued that 
this path leads naturally to non-Gaussian or non¬ 
stationary multivariate random fields, and they fur¬ 
ther made the link with physics-constrained cross¬ 
covariance models. 

Although SLR advocated the SPDE approach as a 
computationally efficient reformulation of univariate 
and multivariate Gaussian random fields, there are 
also limitations that can only be exacerbated in the 
multivariate setting. For example, the smoothness 
parameter of the Matern covariance function is re¬ 
stricted to certain values. Moreover, this methodol¬ 
ogy requires a strong background in numerical anal¬ 
ysis techniques which is not common among statis¬ 
ticians. Thus, we view the multivariate kernel con¬ 
volution approach as a compromise where physics- 
based information can be incorporated, without con¬ 
fronting the numerous difficulties in implementing a 
SPDE framework. Nevertheless, the route of SPDEs 
offers many interesting and challenging avenues for 
future research. 

SLR end their discussion with an important point 
about the inflation of the number of parameters in 
cross-covariance models as the number of variables 
increases. This is indeed a challenging issue when 
applying likelihood methods and various authors 
have resorted to composite likelihood approaches, 
although more investigations in this area are war¬ 
ranted. Finally, the problem of parameter identih- 
ability under infill asymptotics in the multivariate 
setting is mentioned and we see no reason for this 
effect to disappear in this context. 

The discussion of Zhang and Cai (hereafter, ZC) 
centers on trying to understand why and when cok- 
riging does not always outperform kriging. This is 
a very relevant topic given the numerical results 
provided in our two data analyses. To address this 
topic, ZC start by deriving sufficient conditions for 
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the equivalence of Gaussian measures for a bivariate 
Matern cross-covariance model with common length 
scale and smoothness parameters, but allowing for 
different marginal variances. Using ZC’s notation, 
an example of two bivariate Gaussian measures that 
satisfy their sufficient conditions is given by v = 
1; 011,1 = 022,1 = 1, 011,2 = 022,2 = 2, < 712,1 = 1/2, 
012,2 = 1, «i = 2, and a 2 = 1- 

The example provided by ZC where cokriging is 
equivalent to kriging falls into the category of au- 
tokrigeability (Wackernagel, 2003, page 149). A vari¬ 
able is autokrigeable with respect to a set of vari¬ 
ables if the kriging of this variable is equivalent to 
the cokriging. A trivial case is when all variables are 
uncorrelated. Another case, as illustrated by ZC, 
is when the cross-covariance function is separable 
(also called intrinsically correlated in geostatistics). 
Generalizations of the concept of autokrigeability to 
various simplifications of large cokriging systems by 
means of screen effects were investigated by Subra- 
manyam and Pandalai (2008); see also Furrer and 
Genton (2011) for related methods to handle highly 
multivariate spatial data. 

ZC conclude by investigating a situation where the 
auxiliary variable is observed at more locations than 
the predicted variable, leading to the cokriging pre¬ 
dictor being more efficient than the kriging predictor 
as a function of the co-located correlation coefficient, 
r, in a separable exponential cross-covariance func¬ 
tion setting. Interestingly, their formula (9) shows 
that the mean squared prediction error of the cok¬ 
riging predictor can be reduced to at most 1/2 of 
that of the kriging predictor in this particular case. 
We conjecture that the factor 1/2 could be further 
reduced by either adding more observation points in 
O or by adding more variables p 3 (s),..., Y p ( s) in a 
similar cross-covariance function framework. 

Finally, one last issue that remains is the avail¬ 
ability of statistical software for implementing, sim¬ 
ulating and estimating multivariate models. Indeed, 
the choice of a particular model to use in any ap¬ 
plication requires fair expertise in spatial statis¬ 
tics, and we expect these models to become more 
mainstream as well-documented software becomes 
more available. There are some promising pack¬ 
ages available currently that end users should be 
aware of, all of which are available in R. Whereas 
RandomFields contains a large number of multivari¬ 
ate models and can perform high resolution sim¬ 
ulation of these (Schlather et ah, 2014), spBayes 


is geared toward Bayesian analyses of hierarchical 
multivariate models (Finley, Banerjee and Gelfand 
(2015)), and gstat contains some variogram-based 
estimation routines for the linear model of coregion¬ 
alization (Pebesma (2004)). 
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